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' The close agreement between the predictions of dynamical general relativity 

for the radiated power of a compact binary system and the observed orbital 
CD ■ decay of the binary pulsars PSR B1913+16 and PSR B1534+12 allows us to 

bound the graviton mass to be less than 7.6x 10~ 20 eV/c 2 with 90% confidence. 
This bound is the first to be obtained from dynamic as opposed to static field 
relativity. The resulting limit on the graviton mass is within two orders of 
magnitude of that from solar system measurements, and can be expected to 
i . improve with further observations. 
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I. INTRODUCTION 



General relativity assumes that gravitational forces are propagated by a massless gravi- 
ton. Current experimental limits on the graviton mass are based on the behavior of static 
gravitational fields. In particular, a nonzero graviton mass m would cause the gravitational 
potential to tend to the Yukawa form r~ 1 e~ mr , effectively cutting off gravitational interac- 
tions at distances greater than the Compton wavelength mT x of the graviton. The absence 
of these effects in the solar system |1| and in galaxy and cluster dynamics 0Q| thus provides 
an upper limit on m. 

In the dynamical regime, a nonzero graviton mass would produce several interesting 
effects. These include extra degrees of freedom for gravitational waves (e.g., longitudinal 
modes), and propagation at the frequency-dependent speed 

v=y /l- m 2 /oj 2 . (1.1) 

Recently, Will [|j] and Larson and Hiscock || have proposed techniques for examining the 
latter effect with future gravitational- wave interferometer observations to place a limit on m. 
Here we present a new method of bounding the graviton mass, which makes use of existing 
binary pulsar observations. Our technique is based on the agreement between the observed 
orbital decay of the binary pulsars PSR B1913+16 and PSR B1534+12 and the predictions 
of general relativity ||-f§. This is the first bound on m from dynamic-field relativity to be 
accessible with existing observational data, and it provides a limit that is independent of 
the Yukawa bounds. 

The idea is quite simple. Consider the Hulse- Taylor binary pulsar, PSR B1913+16, of 
which the observed decay rate coincides with that expected from relativity to approximately 
0.3%. A nonzero graviton mass would upset this remarkable agreement by altering the 
predicted orbital decay.Q This implies an upper limit on the graviton mass. A crude estimate 
of this bound is quickly obtained from dimensional analysis. For a system with characteristic 
frequency to one expects the effects of a graviton mass to appear at second order in m/u, as 
in ( |1 . 1|) . For gravitational waves at twice the orbital frequency of PSR B1913+16, requiring 
(m/oj) 2 < 0.003 implies an upper limit of order 10 _20 eV/c 2 . This is comparable to the best 
limit from solar system observations, m < 0.44 x 10 -21 eV/c 2 |Q. The purpose of this paper 
is to refine and make rigorous this estimate. 

In section || we discuss linearized general relativity with a massive graviton. The field 
equation and the effective stress tensor for the metric perturbations (gravitational waves) 
are found. In section we solve the field equations using Fourier techniques, and de- 
rive the gravitational-wave luminosity of a general slowly moving periodic source when the 
graviton is massive. We apply this result to the observed orbital decay of the binary pulsars 
PSR B1913+16 and PSR B1534+12 to obtain an upper limit on the mass of the graviton 
in section [TV], and conclude with some brief comments in section M. 



1 Corrections to other characteristics of the system are negligible by comparison on dimensional 
grounds: (mr) 2 = (m / 'uj) 2 [v / c) 2 , (mM) 2 = (m/Lo) 2 (v/c) 6 , where v/c = O(10~ 3 ) for these binary 
systems. 
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II. LINEARIZED GENERAL RELATIVITY WITH A MASSIVE GRAVITON 



In linearized general relativity one writes the metric as a perturbation of the Minkowski 
metric: 



IMU 1 



(2.1) 



We adopt the convention that indices of h^ v are raised and lowered using the Minkowski 
metric; e.g., 



(2.2) 



The linearized theory is defined by substituting (|2.1|) into the Einstein action, expanding 
in powers of h^ u , and keeping only terms up to second order in h^ u (giving field equations 
linear in h^ v ). 

We wish to examine an extension of linearized general relativity which includes a mass 
term for the graviton. We choose the unique mass term for which the wave equation of 
the linearized theory takes the standard form with an /i-independent source, and for which 
the predictions of massless general relativity are recovered by setting m — > at the end of 
the calculations (see |p|,|TQ|| and Appendix 0). Following the procedure described above, we 
arrive at the action 



64vr 



d x 



h^ x h^ x - 2h ll fU lX x + 2h,j v h 



327r/v,T^ + m 2 (V^ - -h 2 ) 



[2.3a) 



where 



h = h\. (2.3b) 

The first five terms are the linearized Einstein action and the stress tensor source for the 
metric perturbations, while the last term is our (phenomenological) choice of mass term [|§ . 
Linearized general relativity is regained by setting m = 0. At linear order the stress tensor 
is assumed to be independent of h^ u and conserved: 

V = 0. (2.4) 

The field equations arise from requiring the action to be invariant under variations of 
the metric perturbation; one finds 

° V - K X M - + Kv» + ^uh p \ pa 

- V^ ah - w 2 (V - -rj^ti) = -16nT^ . (2.5) 

This rather cumbersome equation simplifies considerably when expressed in terms of the 
trace- reversed metric perturbations hn U , defined by 
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h/a, = h,,,, - -jp^h . (2.6) 



The conservation of the stress tensor requires the divergence of both sides of (|2.5|) to vanish. 
This implies that the mass term itself must have vanishing divergence: 

^"V = 0. (2.7) 

This is equivalent to the Lorentz condition of the massless theory. Here, however, it is not a 
gauge choice; rather, it represents the constraints provided by the equations of motion and 
thus eliminates four of the ten independent h^ u . The remaining six components represent 
true degrees of freedom in the massive theory, which consist of the five helicity states of the 
spin-2 field, plus an additional spin-0 component ||TU| . 



Imposing (|2.7|) , the field equation may be simplified to 

(□ - m 2 ) V = -16vrT^ , (2.8) 

which is the familiar form of the wave equation for a massive field. This will be very 
convenient for calculations of gravitational radiation in the massive-graviton theory. As 
described above, our mass term is the unique choice for which the wave equation takes 
this standard form with an /i-independent source, and for which the predictions of massless 
general relativity are recovered by setting m —>■ at the end of the calculations (see PJlO 



and Appendix |A|). The second property is nontrivial since, for example, the pure spin-2 
theory with five degrees of freedom does not satisfy the classical weak-field tests of general 
relativity asm^O [|T~0| 1 2 ] . 



To analyze the energy content of gravitational waves we need an effective stress tensor 
for metric perturbations. Applying Noether's theorem |13| to the Lagrangian of (p.3a|) we 
find 



1 : {ha^h aP ,--lh^h )U ). (2.9) 



32tt x 2 

Here the brackets denote an averaging over at least one period of the gravitational wave. 
Equation ( |2.9|) is identical in form to the usual effective stress tensor for gravitational waves 
with m = 



III. SOLUTIONS 

In linearized general relativity the field equation ( |2.8| ) with m = has the general solution 

n 

V(t,f) = 4 / dV T^-tf-flZ) . ( 3. 1} 
J \x — x | 

For a massive graviton (|3.1| ) is no longer applicable, since the speed of propagation of the 
gravitational waves is frequency- dependent and so the retarded time t— \x — x'\/v{uj) is 
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different for each frequency component of the wave. We evade this difficulty by solving ( |2.8|) 
in frequency space, dealing with each frequency separately. A similar analysis encompassing 
the radiation of general scalar and vector fields can be found in [131 . 
In the frequency domain, the field equation (|2.8|) becomes 



(V 2 + [uj 2 - m 2 }) h^{u\x) 



-167rT Ml/ (o;|f) 



(3.2) 



where the tilde denotes the Fourier transform and V 2 is the 3-space Laplacian. Equation 
( p. 2D is the inhomogeneous Helmholtz equation; the retarded Green function Gr for this 
equation is 



Gr(uj\x; x') 



4tt\x — x' 



(3.3) 



where 



k = sign(u;)v / ^ 2 — 



(3.4) 



for \u\ > m. (The wavenumber k should not be confused with a spatial index.) The retarded 
solution of (13.20 for fixed uj is then 



h^ v {uj\x) = 16n / d x' Gr{ui\x]x h )T^ lv {ui\x') . 



(3.5) 



In order to evaluate (|3.5|) we make use of the slow- motion approximation, uja <C 1, with 
a the characteristic size of the source. With this assumption, and taking the observation 
point far from the source region (r = \x\ 3> |af|), the Green function Gr may be expanded 
for large r. One finds 



4e 



ikr 



h^iulx) 



+ l -{-ik) 2 ' X ' X ' 



d x' T^iuolx' 



1 + (-ik) 



x ■ x 



1 + O | -, uoa) 
r 



(3.6) 



In the m = case one can write the metric perturbations due to a slowly moving source in 
terms of the mass M, dipole moment Dj, and quadrupole moment Ijk of the source, where 



M= ld 3 xT 00 , 



Dj = J d x T 00 xJ , 
Ijk — / d 3 x Too x-* x k 



(3.7a) 
(3.7b) 
(3.7c) 



We can obtain an analogous result in the frequency domain, using the conservation of the 
stress tensor to write the integral over T^ v in ( |3.6| ) in terms of the multipole moments of 
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the source. In the frequency domain the conservation equation ( |2.4| ) for the stress tensor 
becomes 



- iuT m = OaTci 



—iuTm = diTa 



- 00 — u j ± Gj > — tw-i-oi — ^jJ-tj 

Using these relations and the slow-motion approximation, one can show that 



(3.8) 



h j(u>\x) 
h jk (u>\x) 



4e 



ikr 



r 
r 
r 



x-> 



M+—(-ik)D j + n 
r 2r 



2 v -' 7 '-) 2/ ,/.- 



;r 



-{-iu)Dj - —(-ik)(-iw)I jk 



-iu) 2 I jh 



(3.9) 



where M, Dj, Ijk, are respectively the Fourier transforms or Fourier coefficients of the mass, 
dipole moment, and quadrupole moment of the source. Only the quadrupole terms are 
relevant to us; the mass and dipole moments are constant to linear order in h [the energy 
and momentum carried away by the radiation field are 0(h 2 )], hence M and Dj contain 
only zero-frequency components and will not contribute to the radiation. 

The rate of energy loss by the source can be found by integrating the outward 
gravitational- wave flux over a sphere centered on the source: 



dE 



(3.10) 



Let us assume the source is periodic with period P. Then the metric perturbations h^ u (t, x) 
in the time domain are related to their Fourier components h^ v (uj,x) via 



£ 1 



-iuj-nt 



P Jo 



dt hu U (t, x) e 



where 



u n = n 



2tt 



(3.11) 
(3.12) 

(3.13) 



and the tilde now represents a Fourier coefficient. Substituting ( |3.9|) and ( |3.11|) into the 
expression ( j2.9[ ) for the stress tensor of the gravitational waves the luminosity is found to 
be 



L = L GR + J2 



oo 2 4 

m u„ 



n=l 



+ O (m 4 ) 



(3.14a) 



where 
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E< 

n=l 



.OJr. 



2 
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tr/(a;„ 



(3.14b) 



is the usual general-relativistic expression for the radiated power, tr J is the trace of Jj^, 
and we sum over repeated indices. The quantity in the summation of ( 3.14a[ ) is the first 
correction to the general-relativistic expression for the radiated power due to a small nonzero 
graviton mass. Comparison of this correction to the observed orbital decay in binary pulsars 
PSR B1913+16 and PSR B1534+12 will provide us with a bound on m. 



IV. BINARY PULSARS 



The formula 14| ) for the energy-loss rate of a gravitational- wave source when the 



graviton is massive is easily applied to the orbital decay of binary systems. Consider two 
bodies of masses Mi and M 2 , orbiting in the xy plane with coordinates (di cos(6 l ), di sin(9)), 
(—d 2 cos(#), — d 2 sin(#)). Choosing the origin to be at the center of mass, one has 

* = (4 - lb » 

where d is the orbital separation of the binary components, fi is the system's reduced mass, 
and M is its total mass, 

d = d 1 + d 2 , (4.1c) 

» = (4.1d) 
M = M 1 + M 2 . (4.1e) 

Assuming a Keplerian orbit, the motion is described by 

a(l — e 2 ) , 



d6 _ [Ma(l - e 2 )] 1 / 2 
dt ~ d? 



(4.3) 



where a is the semi-major axis and e is the eccentricity of the orbit. The nonzero quadrupole 
moments of this system are 

I xx = /id 2 cos 2 (6>) , 

h y = lyx = ^d 2 cos(6 l ) sin(6>) , 

I yy = fid 2 sin 2 (6) . (4.4) 

The Fourier transform of the quadrupole moment of Keplerian orbits is known For 
n > 
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[UJr. 



2 



2n 



[J n -2(ne) - 2eJ n _i(ne) + 2eJ n+1 (ne) - J n+2 (ne)} 
[l - e 2 )5 [J n _ 2 (ne) - 2J n (ne) + J n+2 (ne)} , 
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-— [J n - 2 {ne) - 2eJ n _i(ne) + -J n (ne) + 2eJ n+1 (ne) - J n+2 (ne)} , (4.5) 
2n n 



where the J n (x) are Bessel functions of the first kind The moments for n < follow from 



(4.6) 



Combining these quadrupole moments with equation ( p.14 ) provides us with an easy 
means to put a limit on the graviton mass. For example, the orbital decay rate of the 
binary pulsar PSR B1913+16 has been measured and found to be slightly in excess of the 
predictions of general relativity 0. Denote by Pb the measured orbital period of the binary 
system, Pb the measured orbital period derivative ascribed to gravitational radiation, and 
Pgr the instantaneous period derivative expected owing to general- relativistic (i.e., zero 
graviton rest mass) orbital decay. Identify the fractional discrepancy between the observed 
and predicted decay rates: 



A = 



P, 



GR 



P 



(4.7) 



GR 



For a slowly decaying Keplerian binary, the instantaneous period derivative is proportional 
to the energy-loss rate; hence, 



Pb-P 



GR 



L-L 



GR 



P 



GR 



^GR 



(4- 



where L is the gravitational- wave luminosity inferred from Pb, and Lgr is the energy- loss 
rate expected from general relativity. This quantity has been measured for PSR B1913+16 
and PSR B1534+12 (see §,§ and Table |). 

Now suppose that A is due at least in part to a nonvanishing graviton mass (rather than 
simply experimental uncertainties). Combining ( |3.14 ) and ( |4.8| ), this implies an upper limit 
to the squared graviton mass of 



24 

m 2 < — F(e) 
5 



2nh 



c 2 P h 



Pb-P 



GR 



P 



(4.9) 



GR 



where P(e) is a function of the eccentricity, 



Fie) 



,5> 6 

1 71=1 



3/jfc(w n )/j I fc (a; n ) 



tr IiuJn 



12 



E 

n=l 



Ij k {u n )I* k {u n ) 



triYov, 



(4.10) 



These sums can be performed using the techniques of pi, giving 
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The function F(e) is plotted in Figure |l|. Note that F(e) is greater than or equal to unity; 
a nonzero graviton mass increases the energy emission of Keplerian binaries, as one would 
expect from adding extra degrees of freedom to the gravitational field. Figure p] contains 
another lesson, as well. Note that, for binaries of fixed period, stronger bounds arise from 
binaries with smaller eccentricity. This dependence is easily understood. Binaries with large 
eccentricities have strong speed variations, as they move from periastron to apiastron. These 
speed variations lead high-eccentricity binaries to produce the bulk of their radiation in ever 
higher harmonics of the orbital frequency ||16|| . The effects of a non-zero graviton mass are 



more pronounced for lower- frequency gravitational waves, as in equation (1.1). As a result, 



the ideal system for bounding the graviton mass is a binary with a large orbital period and 
small eccentricity (a weak emitter of gravitational waves), but which still has a measurable 
inspiral rate. 

Equation ( [4.9| ), relating the squared graviton mass to the fractional discrepancy A in 
the period derivative [or equivalently by ( |4.8| ) the fractional discrepancy in the luminosity], 
assumes that this discrepancy is known exactly. In fact, the period-derivative discrepancy is 
known only up to the errors associated with the measured changes in the binary period and 
the acceleration of the binary relative to Earth. In practice, the one-sigma uncertainty in A 
(which is listed in Table |) is of the same order as the measured discrepancy and must be 
accounted for, as the actual A could reasonably be expected to differ from the value derived 
from the measurements by one or more standard deviations. Consequently, we must describe 
the actual upper-limit on the mass statistically. In the absence of detailed information we 
assume the measured discrepancy A to be normally distributed about its unknown actual 
value [given by the equality in ( |4.9| ) with unknown m 2 ], and with standard deviation as 
given in Table |. In our model we relate the discrepancy to the squared graviton mass, 
which must be non-negative. Referring to Table X], which lists the 90% unified upper 
limit/confidence intervals for the non-negative mean of a univariate normal distribution 
based on a measured sample from the distribution, we calculate the 90% upper limit on the 
(non-negative) graviton mass, which is given in the final row of Table |. 

The best single limit on the graviton mass, m < 6.4 x 10~ 20 eV/c 2 , comes from the 
observations of PSR B1534+12. This is despite the larger uncertainty in the measured 
luminosity discrepancy, compared to PSR B1913+16, because the luminosity discrepancy 
for PSR B 1534+ 12 is negative. A negative discrepancy, taken as exact, would correspond 
to a negative graviton mass, which is unphysical. Correspondingly, a negative measured 
discrepancy is particularly unlikely to arise from a positive m 2 compared to a vanishing m 2 , 
which leads to a tighter upper limit. 

We may combine the two observed discrepancies to find a single upper bound on the 
graviton mass. Each observation k results in a discrepancy A& and an associated one-sigma 
uncertainty in the estimated discrepancy <JA,k- These in turn are related, through equation 
( |4.9| ), with measurements of m 2 kl together with associated one-sigma uncertainties o^. The 
quantity 

m 2 = — — - — - (4.12a) 
1 + P V ; 
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is then a normally distributed random variable whose mean is the squared graviton mass 
and whose variance is 

<x 2 = ^ (4.12b) 



(1+PY 



Choosing 



2 

0"! ' 



P=[f 2 ) (I- 12c-) 

minimizes the variance of m 2 : 

a 2 = 4^L. (4.12d) 
a 2 + a 2 

Referring to Table | and [|I7], Table X], the corresponding limit on the graviton mass from 
the combined observations of PSR B1913+16 and PSR B1534+12 is thus 

m 90 % < 7.6 x KT 20 eV/c 2 . (4.13) 



V. DISCUSSION 

Table | gives the relevant parameters and the corresponding graviton mass bounds for 
the two binary pulsars whose gravitational- wave induced orbital decay has been measured, 
PSR B1913+16 and PSR B1534+12 ||||. The graviton mass bounds from the timing 
observations of each system are very similar, and about two orders of magnitude weaker 
than the Yukawa limit obtained from solar-system observations, m < 4.4 x 10~ 22 eV/c 2 [|TJ. 
Both of these bounds are, in turn, several orders of magnitude weaker than that provided 
by observations of galactic clusters, m < 2 x 10~ 29 eV/c 2 @||, though we regard these 
galactic cluster bounds as less robust, owing to their reliance on assumptions about the 
dark matter content of the clusters, for example. In contrast, the bound obtained here is 
very straightforward and involves few assumptions, making it less prone to error: the chief 
assumption that we have made is the form of the effective mass term for the graviton, which 
- while not unique — is natural. Furthermore, any other mass term would be expected 
from dimensional arguments to yield similar results. 

We have assumed that only measurement errors enter into the determination of the 
intrinsic binary period decay rate P&. In fact, the determination of this rate requires an 
estimate of the acceleration of the binary system, which is principally toward the galactic 
center || . This, in turn, depends on an accurate distance measurement to the binary system, 
which can be difficult to make. A systematic error in this distance estimate leads directly to 
an error in the estimated acceleration of the binary and, in turn, to an error in the P& ascribed 
to gravitational radiation induced decay of the binary system. The large uncertainty in the 
discrepancy A associated with PSR B1534+12 may well be due to an underestimate of the 
distance to this binary system in which case the bound on m 2 would be even tighter. 
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The bound described here arises from the properties of dynamical relativity, making 
it conceptually independent of either the solar system or galactic cluster bounds on the 
graviton mass, which are based on the Yukawa form of the static field in a massive the- 
ory. Furthermore, we expect improvement in the bounds from any given pulsar system 
as observations improve the accuracy of the measured fractional discrepancy in the period 
derivative. For example, when the observations of PSR B1534+12 improve limits on A to 
the same level as is observed today for PSR B1913+16, the corresponding single-system 
bound on the graviton mass should improve to approximately 2 x 10~ 20 eV/c 2 . 

The field of gravitational-wave detection is new. We are only just now learning to exploit 
the opportunities it is creating for us. Within the next year, several large ground-based 
interferometric detectors will begin full operation flq-20i, and existing cryogenic acoustic 



detectors [ZTj -|53] will see significant improvements in sensitivity. Within the next decade we 
should see further enhancements in the capability of all these instruments p5^p7|, and the 
deployment of the space-based interferometric detector LISA (Laser Interferometer Space 
Antenna) |28| , 29] . As gravitational- wave observations mature, we can expect more and 
greater recognition of their utility as probes of the character of relativistic gravity. The 
opening of the new frontier of gravitational-wave phenomenology promises to be an exciting 
and revealing one for the physics of gravity. 
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APPENDIX A: CHOICE OF MASS TERM 



The most general mass term possible for the linearized action ( 2.3a|) is proportional to 

i 

2 



[h fliy h fiu — n(h u u ) 2 ], with k an arbitrary constant. Here we demonstrate that k — \ is the 



unique choice possessing both of the following properties: 

1. the field equations for the metric perturbations can be written in the standard form 

(□-m 2 )V = -16vrT;«, (Al) 



where the source is a local function of the stress tensor and is independent of h^; 
and 

2. taking the limit m — > in the massive theory recovers the predictions of general 
relativity. 

The first property is practical, while the second is necessary for agreement with experiment. 
The field equation for for general n is 
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The divergence of both sides of ( |A2| ) must be equal, implying 



(A3) 



Taking the trace of the field equation and using this divergence condition gives the trace 
condition 



2(1 - K)Uh + (1 - AK)m 2 h = 16ttT 



(A4) 



We see that h can only be written as a local function of the stress tensor if k — 1. 
Substituting the trace and divergence conditions into the field equation gives 



[D - m 2 )h 



2\ ht iv 



-16tt ( --r]^T x x ) + (2k - 1) 



2 ' 



(A5) 



which is of the desired form (|A1|) except for the term in square brackets. The latter can be 
removed only for two special values of k. For k = \ the coefficient vanishes, leaving 



;□ - mr)k 



(A6) 



which is equivalent to (|2.8|) . For k — 1 (the Pauli-Fierz mass term used by Boulware and 
Deser ||10|| ) we can use the trace condition ( |A4| ) to rewrite the term in square brackets as a 
local function of the stress tensor, yielding 



'□ -m 2 W v 



1 



- 1 6 77 ( -^ V T X X 



1 



3m 2 



(A7) 



which is also of the desired form (|A1|) . It is well known, however, that the predictions of 
the k = 1 theory do not reduce to those of general relativity for m — > 0: this is the van 
Dam-Veltman-Zakharov discontinuity |1T| , |I2 |. We are thus led to the choice k = | and the 
massive graviton theory described by ( |2.3a| ). 
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TABLES 



TABLE I. Orbital parameters and corresponding graviton mass bound from the two binary pul- 
sar systems whose gravitational wave induced orbital decay has been measured. Pulsar parameters 
are taken from [ 



; see also [30|. One-sigma uncertainties are quoted for A. 





PSR B1913+16 


PSR B1534+12 


Period 




27907 s 


36352 s 


Eccentricity 




0.61713 


0.27368 


A 




0.32% ± 0.35% 


-12.0% ± 7.8% 


Graviton mass 9 


0% upper bound 


9.5 x 10" 20 eV/c 2 


6.4 x 10~ 20 eV/c 2 



15 




16 



